﻿#########################################################################################
#																						#
# Calculate a variety of outputs from a single wave data file							#
# or from a group of wave data files.													#
# Outputs:																				#
#	Single wave:																		#
#		2D & 3D plots of height, pressure, (experimental) P_x and P_y					#
#		streamlines																		#
#   Batch:																				#
#		same as above for each wave														#
#		an XML database file															#
#		(experimental) pressure RMS around center of the wave (sorted by Q or date)		#
#		amplitude  (sorted by Q or date)												#
#																						#
# Matt Ball (ball@brown.edu)															#
# 2009																					#
#																						#
#########################################################################################

import sys, Gnuplot, os, time
from time import strftime, localtime
from PyTrilinos import Epetra
from numpy import arange, std, array, zeros, ones, average
from scipy.integrate import quad, trapz
from math import pi, tanh

def saveArray(A, path, title, sortType, x_labels):
	"""
	For saving plots of (sortable) "scatterplot" functions, like amplitude & pressure rms
	
	A        -- a 2 x n 2D array of (x,y) values to plot
	path     -- the directory path to save the plot to
	title    -- the title of the plot file, excluding '.png'
	sortType -- either 'Q' or 'date' (not actually used for sorting within this function)
	x_labels -- either None or an array of strings to use as the xtic labels
	"""
	
	g = Gnuplot.Gnuplot()
	# since enhanced mode interprets underscores as a control character
	enhancedTitle = title.replace("_","\\\_")
	g('set title \"' + enhancedTitle + '\"')
	g('set term png size 800,600 enhanced font \"Times New Roman,12\"')
	g.xlabel(sortType)
	g.ylabel(title)
	if not (sortType == None):
		title = title + "_by_" + sortType
	output = '\"' + path + "/" + title + '.png\"'
	output = output.replace("files","pressure")
	g('set output ' + output)
	
	if(x_labels != None):
		xtics_string = ""
		for i in range(len(x_labels)):
			xtics_string = xtics_string + '\"' + x_labels[i] + '\"' + str(i) + ', '
			if (i == len(x_labels)-1):
				# strip off the extra ', '
				xtics_string = xtics_string[:-2]
		
		g('set xtics (' + xtics_string + ')')
		g('set xtics nomirror out rotate')
	
	g.plot(A)
	
def plotStreamlines(A):
	"""
	Takes a height array (2D) and generates a plot of the streamlines for that wave
	"""
	
	title = "streamlines"
	g = Gnuplot.Gnuplot()
	
	numLines = 20 # the number of streamlines shown
	
	rows = A.M()
	cols = A.N()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	#g.xlabel('q')
	#g.ylabel('height')
	g('set term png size 800,600 enhanced font \"Times New Roman,12\"')
	output = '\"' + fileName + '_' + title + '.png\"'
	output = output.replace("files","pressure")
	g('set output ' + output)
	g('set data style dots')
	
	sls = zeros((rows*numLines,2))
	
	col = 0
	for i in range(rows*numLines):
		sls[i][0] = x[i%rows]
		sls[i][1] = A[i%rows][col]
		if (i % rows == 0):
			col += (cols/numLines)
	
	g.plot(sls)
	
def plotMat(A,title):
	"""
	Takes an array (2D) and generates a 2D colormap plot
	saveToDisk is a global boolean - 'False' for single wave, 'True' for batch 
	"""
	# 'title' is "height", "pressure", etc.
	g = Gnuplot.Gnuplot()
	rows = A.M()
	cols = A.N()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	g('set pm3d map')
	g.xlabel('q')
	g.ylabel('p')
	
	if (saveToDisk):
		g('set term png size 800,600 enhanced font \"Times New Roman,12\"')
		output = '\"' + fileName + '_' + title + '_2D.png\"'
		output = output.replace("files","pressure")
		g('set output ' + output)
	
	g.splot(Gnuplot.GridData(A,x,y, binary=0))
	
def plotMatForAnimation(A, title, time):
	"""
	Takes an array (2D) and generates a 2D colormap plot named suitably for
	creating an animated gif
	"""
	# 'title' is "height", "pressure", etc.
	g = Gnuplot.Gnuplot()
	rows = A.M()
	cols = A.N()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)*-1
	g('set xrange [-pi:pi]')
	g('set pm3d map')
	g('unset colorbox')
	g('set format y \"-%g\"')
	
	if (saveToDisk):
		g('set term png size 1024,768 enhanced font \"Times New Roman,12\"')
		output = '\"' + fileName + '_' + title + '_2D.png\"'
		output = output.replace("files","animation")
		output = output.replace("sol.",str(int(time)) + "_sol.")
		g('set output ' + output)
	
	g.splot(Gnuplot.GridData(A,x,y, binary=0))
	
def plotMat3D(A,title):
	"""
	Takes an array (2D) and generates a 3D surface plot
	saveToDisk is a global boolean - 'False' for single wave, 'True' for batch 
	"""
	# 'title' is "height", "pressure", etc.
	rows = A.M()
	cols = A.N()
	g = Gnuplot.Gnuplot()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	g('set pm3d')
	g.xlabel('q')
	g.ylabel('p')
	g('unset surface')
	g('set view ,225,1,1')
	
	if (saveToDisk):
		g('set term png size 800,600 enhanced font \"Times New Roman,12\"')
		output = '\"' + fileName + '_' + title + '_3D.png\"'
		output = output.replace("files","pressure")
		g('set output ' + output)
	
	g.splot(Gnuplot.GridData(A,x,y, binary=0))
	
def makePressureDirs(path):
	"""
	Walk a directory tree and make a directory called 'pressure'
	whenever there is a directory called 'files', ignore if already present 
	"""
	for dirpath, dirnames, filenames in os.walk(path):
		for dir in dirnames:
			if (dir == "files"):
				path = os.path.join(dirpath,"pressure")
				try:
					os.mkdir(path)
				except (OSError):
					pass
	
def readFile(file, name, saveIt):
	"""
	The main function which operates on a wave. It reads a wave data file line-by-line
	and generates a variety of plots. To control which plots are generated, (un)comment
	the appropriate lines under the '#height', '#pressure', '#streamlines', etc. comments.
	"""
	
	global p0
	global fileName
	global saveToDisk
	global Q
	global gamma
	global type
	global A # used for tanh
	global b
	
	fileName = name
	saveToDisk = saveIt
	
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	# the size of gamma here is *rows* since the matrix gets transposed after being read in
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
		if (b == 0):
			type = "exp"
			gamma = A*ones(rows) # only for b=0 since gamma = A * exp[b*t] in general (for 'exp')
		else:
			# VERY SPECIAL CASE: tanh - NOT EXP at all!
			type = "tanh"
			gamma = zeros(rows) # to get delta_p in myGamma(...)
			
	elif(type == "step\n"):
		p1 = float(file.readline())
		p2 = float(file.readline())
		a0 = float(file.readline())
		a1 = float(file.readline())
		a2 = float(file.readline())
		type = "step"
		gamma = step_gamma(rows, p0, p1, p2, a0, a1, a2)
	
	file.readline() # blank line

	h = Epetra.SerialDenseMatrix(cols,rows,False)
	for i in range(rows):
		for j in range(cols):
			try:
				h[j][i] = float(file.readline())
			except ValueError:
				print "ValueError occurred somewhere near (probably after) line", `(i*cols)+j`
				print "The file is probably shorter than expected, sorry :("
				return
	
	h = centerWave(h)

	# height
	plotMat3D(h, "height")
	plotMat(h, "height")
	#plotMatForAnimation(h, "height", os.stat(fileName).st_mtime)
	
	# pressure
	pressure = P(h)
	plotMat3D(pressure, "pressure")
	plotMat(pressure, "pressure")
	
	# streamlines
	plotStreamlines(h)
	
	# P_x and P_y
	#pressure_x = P_x(h, pressure)
	#pressure_y = P_y(h, pressure)
	#plotMat3D(pressure_x,"P_x")
	#plotMat3D(pressure_y,"P_y")
	#plotMat(p_x,"P_x")
	#plotMat(p_y,"P_y")
	
def step_gamma(numElts, p0, p1, p2, a0, a1, a2):
	"""
	Build a 1D array of \Gamma(p) for step
	Should look like [a2, a2, a2, ..., a1, a1, a1, ..., a0, a0, a0]
	"""
	delta_p = abs(p0)/numElts
	temp = zeros(numElts)

	for i in range(int(abs(p2)/delta_p)):
		temp[i]= a2
	for i in range(int(abs(p2)/delta_p), int(abs(p1)/delta_p)):
		temp[i] = a1
	for i in range(int(abs(p1)/delta_p), int(abs(p0)/delta_p)):
		temp[i]= a0
	
	return temp

def myGamma(numElts,delta_p):
	"""
	Calculate the \Gamma array from the variable named gamma in readFile
	This is named myGamma since Python already has a gamma function defined
	"""
	Gamma = zeros(numElts)
	
	if (type == "exp") or (type == "step"):
		for i in range(numElts):
			Gamma[i] = trapz(gamma[:i+1], None, delta_p)
	
	elif (type == "tanh"):
		
		for i in range(numElts):
			# ExpWaveProb.h - line 29 for the formula (their Gamma is my ɣ)
			# 		// double Gamma(double p) { return A*tanh(B*(-abs(p)+1)); }
			Gamma[i] = quad(lambda p: A*tanh(b*(-abs(p)+1)), 0, delta_p*(i+1))[0]
	
	return Gamma
	
def centerWave(h):
	"""
	Centers the wave in the matrix
	"""
	# the row at p = 0, iterate over q
	rows = h.M()
	cols = h.N()
	
	indexOfCenter = 0
	max = 0
	for i in range(rows):
		if (h[i][0] > max):
			max = h[i][0]
			indexOfCenter = i
	
	shiftBy = (rows/2) - indexOfCenter # this must be an int
	
	newWave = Epetra.SerialDenseMatrix(rows,cols,False)
	
	for i in range(rows):
		for j in range(cols):
			newWave[i][j] = h[(i-shiftBy)%rows][j]
	
	return newWave

def dq(h):
	"""
	Returns the partial derivative w/r/t q of whatever function h is
	Note, this uses periodic boundary condition to evaluate centered difference at edges
	"""
	rows = h.M()
	cols = h.N()
	
	delta_q = 2*pi/rows
	
	dq = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			# don't need to worry about negative indices, python takes care of them
			dq[i][j] = h[(i+1)%rows][j] - h[i-1][j]
			
	return dq/(2*delta_q)

def dp(h):
	"""
	Returns the partial derivative w/r/t p of whatever function h is
	Note, this uses uses forward or backward difference at edges
	"""
	rows = h.M()
	cols = h.N()
	
	delta_p = abs(p0)/cols
	
	dp = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			if (j==0):
				# forward
				dp[i][j] = (h[i][j+1] - h[i][j])/delta_p
				
			elif (j==cols-1):
				# backward
				dp[i][j] = (h[i][j] - h[i][j-1])/delta_p
				
			else:
				# general case
				dp[i][j] = (h[i][j+1] - h[i][j-1])/(2*delta_p)
			
	return dp
	
def P(h):
	"""
	Returns the pressure, given the height
	"""
	rows = h.M()
	cols = h.N()
	
	hq = dq(h)
	hp = dp(h)
	
	#g = 9.806 # For comparison to legacy code ONLY
	g = 9.816
	
	delta_p = abs(p0)/cols
	Gamma = myGamma(cols, delta_p)
	
	P = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			P[i][j] = - .5*(1 + hq[i][j]**2)/(hp[i][j]**2) - g*h[i][j] - Gamma[j] + .5*Q
	
	return P
	
def P_x(h,P):
	"""
	Returns \partial pressure / \partial x, given height and pressure
	Takes P in addition to h to avoid recalculating too much
	"""
	rows = h.M()
	cols = h.N()
	
	hq = dq(h)
	hp = dp(h)
	
	Pq = dq(P)
	Pp = dp(P)
	
	Px = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			Px[i][j] = Pq[i][j] - (Pp[i][j])*(hq[i][j])/(hp[i][j])
	
	return Px

def P_y(h,P):
	"""
	Returns \partial pressure / \partial y, given height and pressure
	Takes P in addition to h to avoid recalculating too much
	"""
	rows = h.M()
	cols = h.N()
	
	hp = dp(h)
	
	Pp = dp(P)
	
	Py = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			Py[i][j] = Pp[i][j]/hp[i][j]
	
	return Py
	
def relpath(path, start='.'):
	"""
	Methods borrowed from Python v2.6 (kind of hacked) to make some other stuff work
	"""
	
	# Return the longest prefix of all list elements.
	def commonprefix(m):
		"Given a list of pathnames, returns the longest common leading component"
		if not m: return ''
		s1 = min(m)
		s2 = max(m)
		for i, c in enumerate(s1):
			if c != s2[i]:
				return s1[:i]
		return s1
	
	"""Return a relative version of a path"""

	if not path:
		raise ValueError("no path specified")

	sep = '/'
	pardir = '..'
	
	start_list = os.path.abspath(start).split(sep)
	path_list = os.path.abspath(path).split(sep)

	# Work out how much of the filepath is shared by start and path.
	i = len(commonprefix([start_list, path_list]))

	rel_list = [pardir] * (len(start_list)-i) + path_list[i:]
	if not rel_list:
		return curdir
	return os.path.join(*rel_list)

def getPressureStdev(file, percentAroundCenter):
	"""
	Returns the standard deviation of the pressure along p=0
	This is function is very inefficient since it computes the pressure over the entire wave! (optimize later)
	percentAroundCenter in [0,100] changes the size of the neighborhood around the center
	"""
	
	global p0
	global type
	global gamma
	global type
	global Q
	
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
		type = "exp"
		gamma = A*ones(cols) # only for b=0 since gamma = A * exp[b*t] in general (for 'exp')
	elif(type == "step\n"):
		# confirm which is which
		p1 = float(file.readline())
		p2 = float(file.readline())
		a0 = float(file.readline())
		a1 = float(file.readline())
		a2 = float(file.readline())
		type = "step"
		gamma = step_gamma(cols, p0, p1, p2, a0, a1, a2)
	
	file.readline() # blank line
	
	temp = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			temp[i][j] = float(file.readline())

	# h
	h = Epetra.SerialDenseMatrix(cols,rows,False)
	for i in range(cols):
		for j in range(rows):
			h[i][j] = temp[j][i]

	h = centerWave(h)
	pressure = P(h)

	lowerBound = int(rows/2 - (rows*percentAroundCenter)/200)
	upperBound = int(rows/2 + (rows*percentAroundCenter)/200)
	numElts = upperBound - lowerBound

	topRow = zeros(numElts)
	for i in range(numElts):
		topRow[i] = pressure[i+lowerBound][0]
	
	return std(topRow)
	
def getAmplitude(file):
	"""
	Calculates the amplitude of a single wave
	"""
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
	elif(type == "step\n"):
		# p1, p2, a1, a2, a3
		file.readline()
		file.readline()
		file.readline()
		file.readline()
		file.readline()
	
	file.readline() # blank line
	
	# cols, not rows, because of the transpose issue
	topRow = zeros(cols)
	for i in range(cols):
		topRow[i] = file.readline()
	
	return max(topRow) - average(topRow)
	
def getQ(file):
	"""
	Grabs Q from a single wave
	"""
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	return Q

def byDate():
	"""
	Outputs pressure RMS and amplitude plots, sorted by date
	"""
	pressureRmsAndAmplitude("date")
	
def byQ():
	"""
	Outputs pressure RMS and amplitude plots, sorted by Q
	"""
	pressureRmsAndAmplitude("Q")

# amplitude and pressure RMS as functions of sortType (either date or Q)
def pressureRmsAndAmplitude(sortType):
	"""
	Outputs pressure RMS and amplitude plots
	"""
	
	if (sortType != "date") and (sortType != "Q"):
		# in case the specified sort type isn't one of the implemented ones
		# check once at the beginning so I don't have to raise an exception later
		raise Exception, "sort type '" + sortType + "' not available"
	
	if ((sys.argv[2][-6:] == "files/") or (sys.argv[2][-5:] == "files")):
		print "Input directory must be a directory containing 'files' and not 'files' itself."
	else:
		makePressureDirs(sys.argv[2])
		numPoints = 0
	
		# count the number of data points we'll be generating
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
			for f in filenames:
				if (dirpath[-6:] == "/files"):
					if (f[:4] == "sol."):
						numPoints = numPoints + 1
		
		print "Computing pressure RMS and amplitude values for", numPoints, "files (sorted by", sortType + ")"
		rms = zeros((numPoints,2))
		amp = zeros((numPoints,2))
		
		x_array = zeros(numPoints)
		x_labels = None
		
		if (sortType == "date"):
			x_labels = zeros(numPoints, dtype='|S20') # 20 = (max) length of the strings in the array
		
		i = 0
		
		# now actually do it
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
			
			if (sortType == "date"):
				# this only works for date modified, not for date created
				filenames.sort(key = lambda f:os.stat(os.path.join(dirpath,f)).st_mtime)
			#else (sortType == "Q"):
				# filenames is already in alphabetical order, so do nothing
			#	pass
			
			for f in filenames:
				
				if (dirpath[-6:] == "/files"):
					if (f[:4] == "sol."):
						name = os.path.join(dirpath, f)
						
						print "(" + str(i+1) + ")\t" + name
						
						file = open(name,"r")
						
						if(len(sys.argv) == 4):
							rms[i][1] = getPressureStdev(file, int(sys.argv[3]))
						else:
							rms[i][1] = getPressureStdev(file, 100)
						
						file.close()
						file = open(name,"r")
						amp[i][1] = getAmplitude(file)
						file.close()
						
						if(sortType == "Q"):
							file = open(name,"r")
							x_array[i] = getQ(file)
							file.close()
						else:
							# sortType equals "date"
							x_array[i] = i
							x_labels[i] = str(strftime("%Y-%m-%d %H:%M:%S",localtime(os.stat(os.path.join(dirpath,f)).st_mtime)))
							
						i += 1
						
		for j in range(numPoints):
			rms[j][0] = x_array[j]
			amp[j][0] = x_array[j]
		
		pressureTitle = ""
		if (len(sys.argv) == 4):
			pressureTitle = "pressure_rms_" + sys.argv[3] + "%"
		else:
			pressureTitle = "pressure_rms_100%"
		
		saveArray(rms, dirpath, pressureTitle, sortType, x_labels)
		saveArray(amp, dirpath, "amplitude", sortType, x_labels)

# creates/updates the XML database (<waveset>...</waveset>) with entries of the form:
#	<wave class="exp" p0="-2" Q="21.5" amplitude="3.4" path="foo/bar/baz" A="5"
#		resolution_q="500" resolution_p="250" >
#		<graph dims="2" type="pressure" path="foo/bar/graph" />
#	</wave>
# important: path names are relative to the directory this script was called from
# only operates on 'sol.' files found in directories named 'files'
def makeDatabase():
	"""
	Create/update the XML database (<waveset>...</waveset>) named 'database.xml'
	Entries are of the form:
		<wave class="exp" p0="-2" Q="21.5" amplitude="3.4" path="foo/bar/baz" A="5"
			resolution_q="500" resolution_p="250" >
			<graph dims="2" type="pressure" path="foo/bar/graph" />
		</wave>
	Important: path names are relative to the directory this script was called from
	This only operates on 'sol.' files found in directories named 'files'
	"""
	
	print "\nMake sure to run this script from the wave data *root directory* when generating the database.\n"
	
	dbName = "database.xml"
	dbPath = os.path.join(sys.argv[2], dbName)
	
	file = None
	
	if (os.path.exists(dbPath)):
		# there already exists a file called database.xml in the directory spec'd by argv[2]
		YorN = raw_input("Overwrite old database file? ")
		yorn = YorN.lower()
		if (yorn not in ["y", "yes"]):
			print "exiting..."
			return
		else:
			pass
			
	# check that argv[2] is a directory
	if (os.path.isdir(sys.argv[2])):
		# it is
		file = open(dbPath, "w") # w mode overwrites the old file, w+ appends
		
		file.write("<waveset>\n")
		numEntries = 0
		numFailedEntries = 0
		
		# Write the entries
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):

			for f in filenames:
						
				if (dirpath[-6:] == "/files"):
				
					if (f[:4] == "sol."):
						name = os.path.join(dirpath, f)
						
						thisEntry = dbEntry(name)
						# don't write empty entries
						if not (thisEntry == "\t<wave \n\t</wave>\n"):
							file.write(thisEntry)
							numEntries += 1
						else:
							numFailedEntries += 1
			
		file.write("</waveset>")
		file.close()
		print "database.xml is complete, with", numEntries, "waves successfully entered and", numFailedEntries, "waves failed"
		print "Currently, this file must be broken up into db_exp.xml, db_step.xml, and db_tanh.xml for use by the web app."
		return
		
	elif (os.path.isfile(sys.argv[2])):
		# it's a file, so we'll treat it as a single wave to operate on
		print dbEntry(sys.argv[2])
		return
	else:
		print "Invalid path, exiting..."
		return

def dbEntry(filePath):
	"""
	Returns a string representing a database entry for a single file
	"""
	file = open(filePath,"r")
	amp = getAmplitude(file)
	file.close()
	
	entry = "\t<wave "

	try:
		file = open(filePath)
		time = int(os.stat(filePath).st_mtime)
		type = file.readline() # exp or step
		rows = int(file.readline())
		cols = int(file.readline())
		p0 = float(file.readline())
		Q = float(file.readline())
		
		# the size of gamma here is *rows* since the matrix gets transposed after being read in
		if (type == "exp\n"):
			A = float(file.readline())
			b = float(file.readline())
			if (b == 0):
				type = "exp"
				entry += "class=\"" + type + "\" p0=\"" + str(p0) + "\" Q=\"" + str(Q)
				entry += "\" amplitude=\"" + str(amp) + "\" path=\"" + relpath(filePath)
				entry += "\" A=\"" + str(A) + "\" resolution_q=\"" + str(rows)
				entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
				
			else:
				# VERY SPECIAL CASE: tanh - NOT EXP at all!
				type = "tanh"
				entry += "class=\"" + type + "\" p0=\"" + str(p0) + "\" Q=\"" + str(Q)
				entry += "\" amplitude=\"" + str(amp) + "\" path=\"" + relpath(filePath)
				entry += "\" A=\"" + str(A) + "\" B=\"" + str(b) + "\" resolution_q=\"" + str(rows)
				entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
				
		elif(type == "step\n"):
			p1 = float(file.readline())
			p2 = float(file.readline())
			a0 = float(file.readline())
			a1 = float(file.readline())
			a2 = float(file.readline())
			type = "step"
			entry += "class=\"" + type + "\" Q=\"" + str(Q) + "\" p0=\"" + str(p0)
			entry += "\" p1=\"" + str(p1) + "\" p2=\"" + str(p2) + "\" a0=\"" + str(a0)
			entry += "\" a1=\"" + str(a1) + "\" a2=\"" + str(a2) + "\" amplitude=\"" + str(amp)
			entry += "\" path=\"" + relpath(filePath) + "\" resolution_q=\"" + str(rows)
			entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
		
		file.close()
	
	except:
		pass
	
	entry += "\n\t</wave>\n"
	return entry

def batch():
	"""
	Does readFile() on each 'sol.' file found in a directory named 'files' and saves plots to disk
	Prints out a percent completed to the terminal (since batch can take a long time) but it's hacked
	"""
	
	# count the number of data points we'll be generating
	numPoints = 0
	for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
		for f in filenames:
			if (dirpath[-6:] == "/files"):
				if (f[:4] == "sol."):
					numPoints = numPoints + 1
		
	print "Generating plots for", numPoints, "waves. This could take a while."
	YorN = raw_input("Proceed? ")
	yorn = YorN.lower()
	if (yorn not in ["y", "yes"]):
		print "exiting..."
		return
	else:
		pass

	makePressureDirs(sys.argv[2])
	
	# ugly hack to periodically print out % complete. add more elements to this array for more frequent messages.
	hacked = array([True, True, True, True, True, True, True, True, True, True])
	percentComplete = 0
	
	for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
	
		for f in filenames:
					
			if (dirpath[-6:] == "/files"):
			
				if (f[:4] == "sol."):
					try:
						percentComplete += (1./numPoints)*100
						
						for i in range(len(hacked)):
							if (percentComplete >= (i*(100/len(hacked)))) and hacked[i]:
								print str(int(percentComplete))+"% complete"
								hacked[i] = False
								
						name = os.path.join(dirpath, f)
						print name
						file = open(name,"r")
						readFile(file, name, True)
						file.close()
					except:
						pass
							
def single():
	"""
	Runs readFile() on the file passed in to the script
	"""
	name = sys.argv[2]
	file = open(name,"r")
	readFile(file, name, False)
	file.close()
	
def main():

	funcmap = dict(single=single, batch=batch, pRMS_amp_byQ=byQ, pRMS_amp_byDate=byDate, db=makeDatabase)
	def usage():
		print 'Usage: python %s [%s]' % (sys.argv[0], ' | '.join(funcmap.keys())) + ' path [RMS % around center]'
	funckey = (len(sys.argv) >= 3) and (len(sys.argv) <= 4) and (sys.argv[1]) or (None)
	func = funckey in funcmap and funcmap[funckey] or usage	
	func()

if __name__ == '__main__':
    main()